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Free energies of crystalline solids: a lattice-switch Monte Carlo method 

A.D. Bruce, N.B. Wilding & G.J. Ackland 

Department of Physics and Astronomy, The University of Edinburgh 
Edinburgh, EH9 3JZ, Scotland, United Kingdom 

We present a method for the direct evaluation of the difference between the free energies of two 
crystalline structures, of different symmetry. The method rests on a Monte Carlo procedure which 
allows one to sample along a path, through atomic-displacement-space, leading from one structure 
to the other by way of an intervening transformation that switches one set of lattice vectors for 
another. The configurations of both structures can thus be sampled within a single Monte Carlo 
process, and the difference between their free energies evaluated directly from the ratio of the 
measured probabilities of each. The method is used to determine the difference between the free 
energies of the fee and hep crystalline phases of a system of hard spheres. 
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One of the fundamental tasks of theoretical condensed 
matter physics is to understand the observed structures 
of crystalline materials in terms of microscopic models of 
the atomic interactions. The principles involved are well 
known: one needs to evaluate which of the candidate 
structures has the lowest free energy for given (model 
and thermodynamic) parameters. In practice the task 
is rather less straightforward. Conventional Boltzmann 
importance sampling Monte Carlo (MC) methods do not 
yield the free energy Q. It is therefore customary to re- 
sort to integration methods (IM) which determine free 
energies by integrating free-energy derivatives measured 
at intervals along a parameter-space-path connecting the 
system of interest to a reference system whose free energy 
is already known. This procedure has been used widely, 
and with ingenuity |^]. Nevertheless it leaves much to 
be desired. In particular, to determine the difference 
between the free energies of two phases one has to re- 
late each of them separately to some reference system, 
with uncertainties which are not always transparent, and 
which can be significant on the scale of the free energy 
difference of interest. Clearly, one would prefer a method 
which focuses more directly on this difference. The ele- 
ments of such a strategy are to be found in the extended- 
sampling techniques pioneered by Torrie and Valleau j| , 
and recently revitalized in the multicanonical method of 
Berg and Neuhaus 0. The key concept underlying this 
method is that of a configuration-space-path comprising 
the macrostates of some chosen macroscopic property M.. 
The method utilizes a sampling distribution customized 
to even out the probabilities of different .M-macrostates. 
In principle, it allows one to sample along a path (whose 
canonical probability is generally extremely small) cho- 
sen to connect the distinct regions of configuration space 
associated with two phases; the difference between the 
free energies of the two phases can then be obtained di- 
rectly from the ratio of the probabilities with which the 
system is found in each of the two regions. This idea has 
been applied in the investigation of the phase behavior 
of ferromagnets 0] , fluids H , and lattice gauge theories 
||. However its application to structural phase behav- 
ior faces a distinctive problem: finding a path that links 
the regions of configuration space associated with two 
different crystal structures pi] , without traversing regions 
of non- crystalline order, which present problems |q] for 
even multicanonical Monte-Carlo studies. We show here 
how one may construct such a path, and use it for direct 
high-precision measurement of free-energy differences of 
crystal structures. 

The idea is simple; we describe it first in general and 
qualitative terms. The atomic position coordinates are 
written, in the traditions of lattice dynamics, as the sum 
of a lattice vector ||], and a displacement vector. The 
configurations associated with a particular structure are 
explored by MC sampling of the displacements. Given 
any configuration of one structure one may identify a 
configuration of the other, by switching one set of lat- 
tice vectors for the other, while keeping the displacement 



vectors fixed. Such lattice switches can be incorporated 
into the MC procedure by regarding the lattice type as a 
stochastic variable. Lattice switches have an intrinsically 
low acceptance probability, since typically they entail a 
large energy cost. But the multicanonical method can 
be used to draw the system along a path comprising the 
macrostates of this 'energy cost', and thence into a re- 
gion of displacement-space in which the 'energy cost' is 
low, and the lattice switch can be implemented. The net 
result is a MC procedure which visits both structures in 
the course of a single simulation, while never moving out 
of the space of crystalline configurations. The method is, 
we believe, potentially very general. We illustrate it here 
by using it to determine the difference between the free 
energies of the two close-packed structures {fee and hep) 
of a system of hard spheres. This problem has a long 
history JhJ . The difference between the free energies (ef- 
fectively, the entropies) is extremely small, and 

recent IM studies have disagreed on its value [Ojl^ . 
It thus provides an exacting and topical testing ground 
for the method proposed here |L3| . 

We consider a system of N particles with spatial coor- 
dinates {r}. The particles are confined in a fixed volume 
V, with periodic boundary conditions |l4j]. We make the 
decomposition 

fi= Ri+ Ui (1) 

where the vectors Ri,i = 1 . . .N = {R} a define the sites 
of a lattice of type a (here, either fee or hep). Clearly 
there are many transformations that will map one set of 
vectors into the other; the mapping we have chosen is 
explained in Fig. l(a),(b): it exploits the fact that the 
two structures differ only in the stacking pattern of the 
close-packed planes. 

We define a partition function(and free energy) asso- 
ciated with the structure a by |ll| 

Z(N,V,T,a)= I JJ[^]exp [-$({«}, a)] 

J{u}£a i 

= exp[-F a (N,V,T)/kT] (2) 

where $ represents the dimensionless configurational en- 
ergy. In the present context 

*m) = {I ^^'^ w 

where a is the hard sphere diameter. The a-label at- 
tached to the integral in Eq. (|^) signifies that it must in- 
clude only contributions from configurations within the 
subspace associated with the structure a |f(| . 

Consider now the canonical ensemble with probability 
distribution 

P({u}, a | N, V, T) = eXp[ 7^ ({ " } ' a)] (4) 
a J. i , » > z(N,V,T) y ' 
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where Z(N,V,T) = J2 a z {N,V,T,a). The probability 
that the system will be found to have structure a provides 
a measure of the associated partition function: 

P(a \N,V,T)= f Y[[dui]P({u}, a\N,V, T) 

J{u}ea i 

Z(N, V, T, a) 



P(M | N, V, T) oc P(M | N, V, T, {v})e 



-V(M) 



(10) 



Z(N,V,T) 



(5) 



The difference between the free-energies of the two 
structures may be thus be expressed as 

F hcp( N > V > T ) ~ F fcc( N > V > T ) = NkTAf ^kTlnU (6) 



where, 



K 



Z(N,V,T,fcc) = P(fcc\N,V,T) 
Z(N, V, T, hep) P{hcp | TV, V, T) 



(7) 



This identification is useful only if one can devise a 
MC procedure that will actually visit the configurations 
{u}, a with the probabilities prescribed by Eq. (Q). To do 
so one must deal with the ergodic block against lattice 
switches ('updates' of the lattice label, a): almost in- 
variably such a switch maps an accessible configuration 
of one structure onto an inaccessible configuration of the 
other (one which violates the hard-sphere constraint im- 
plied by Eq. (||)). Fig. 1(b) provides an example. The 
resolution is to 6ms the sampling procedure so as to favor 
the occurrence of configurations which transform without 
violating this constraint. To do so we define an overlap 
order parameter 



M({u}) = M({u}, hep) - M({u},fcc) 



(8) 



where M{{u},a) counts the number of pairs of over- 
lapping spheres associated with the configuration {u},a 
(again, see Fig. 1(b)). Since M({u},a) will necessarily 
be zero for any set of displacements {u} actually visited 
when the system has lattice a, the order parameter A4 
is necessarily > (< 0) for realizable configurations of 
the fee (hep) structure. The displacement configurations 
with A4 = are accessible in both structures and thus 
offer no barrier against lattice switches. Accordingly the 
set of ./Vf-macrostates provides us with the required 'path' 
connecting the two phases, through a lattice-switch at 
A4 — 0. To pick out this path we must sample from the 
biased configuration distribution 

P({u}, a | N, V, T, {r)}) cx P({u}, a\N,V, T)e^ M ^^ 

(9) 

where {rj} = r](M),M. = 0, ±1,±2... define a set of 
multicanonical weights Q , which have to be determined 
such that configurations of all relevant .A/f-values are sam- 
pled. Once this is done, one can measure the weighted 
distribution of .M-values, and reweight (unfold the bias) 
to determine the true canonical form of this distribution: 



Finally, the difference between the free energies of the two 
structures may be read off from this distribution through 
the identification (cf Eqs. @,(@)) A/ = TV -1 haft, with 



n 



Em>qP(M\N,V,T) 
J2 m<0 P(M\N,V,T) 



(11) 



We have implemented this procedure to study systems 
of N=216, 1728 and 5832 hard spheres (forming, respec- 
tively 6, 12 or 18 close-packed layers). The volume V 
was chosen such that the fraction of space filled, p, satis- 
fies p/pep = 0.7778 JTJ], where Pep = 0.7404 is the space 
filling fraction in the closest packing limit. The MC pro- 
cedure entails sampling the displacement variables {u} 
and the lattice label a. The variables {u} were updated 
by drawing new values from a top-hat distribution |l§| ] 
and accepting them provided they satisfy the hard sphere 
constraint; the lattice switches were attempted (and ac- 
cepted with probability 1/2) only when the system is in 
the M. = macrostate. The weights (which enable the 
system to reach this special macrostate) were obtained 
using methods explained elsewhere (?J . We allowed typi- 
cally 2 x 10 4 sweeps for equilibration and up to 5 x 10 7 
sweeps for final sampling runs on the largest system size. 
The simulations were conducted on DEC ALPHA work- 
stations using overall some 800 hours CPU time. 

Fig. 2 shows the measured overlap distribution for the 
system of N = 1728 spheres; the inset shows the prob- 
ability on a logarithmic scale, exposing the enormity of 
the entropic 'barrier' (probability 'trough') that the mul- 
ticanonical weighting enables us to negotiate. The dif- 
ference between the free energies of the two structures is 
identifiable immediately and transparently from the ratio 
of the integrated weights of the two essentially gaussian 
peaks. Our results for this system and other system sizes 
are gathered together in Table 1, along with the results 
of other authors. 

From Table 1 it is apparent that the present work 
greatly refines the largely inconclusive results of the origi- 
nal IM study Q . Our results are consistent with -though 
substantially more precise than- very recent IM studies 
of Bolhuis and Frenkel |l^] . They are inconsistent with 
the result reported by Woodcock given that A/ is 
believed to decrease as the density is reduced, towards 
melting |L9]]. While we have not attempted an explicit 
analysis of the finite-size behavior, the close agreement 
between our results for N = 1728 and N — 5832 indi- 
cates that the latter should provide an extremely good 
estimate of the thermodynamic limit, confirming the sta- 
bility of the fee structure at this density. 

Our principal concern here, however, is with the gen- 
eral lessons that can be learned about the method intro- 
duced in this work. The precision we have achieved with 
this method is self-evidently a significant advance on that 
of IM studies. Admittedly, this level of precision has en- 
tailed substantial processing time, principally because of 
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the relative slowness of the diffusive exploration of the 
multicanonically weighted configuration space. But the 
point is that the procedure is practicable [^0| , with a com- 
putational strategy that is, we suggest, less complex and 
more transparent than that of IM. Thus, for example, the 
method described in |g| involves integration (of a mean- 
square displacement) along a parameter-space-path con- 
necting each structure to a reference system, comprising 
an Einstein-model of the same structure; the MC integral 
then has to be combined with the known free energy of 
the Einstein model, a virial correction, and a correction 
to the virial correction, before taking the difference be- 
tween the results for the two structures. The uncertain- 
ties in all the contributions have to be assessed separately. 
By contrast, the present method focuses directly on the 
quantity of interest (the relative weights of the peaks in 
Fig. 2); and the precision with which it is prescribed is 
defined by standard MC sampling theory. 

Finally we comment briefly on the more general ap- 
plicability of the method. For systems other than hard 
spheres, the role of the overlap order parameter is played 
by the energy barrier encountered in the lattice switch; 
the generalization of the weighting procedure should be 
straightforward. It seems unlikely that many problems 
will require the level of precision needed here, where the 
two phases are so finely balanced. However some cir- 
cumstances will not generally prove as favorable. In the 
present case one can readily identify a lattice-to-latticc 
mapping which guarantees no overlaps (high-energy-cost 
interactions) amongst subsets of the atoms (those lying 
within the same close-packed plane). The optimal form 
of mapping may not always be so evident. It may also 
prove advantageous to relax the constraint imposed in 
Eq. (Q) that the coordinates {u} represent identical dis- 
placement patterns in the two structures. At the cost of 
only little extra computational complexity the true dis- 
placements may be represented as structure-dependent 
(even site-dependent) functions of a (still common) set 
of coordinates {u}. It should then be possible to ensure 
that a typical displacement pattern in one structure maps 
onto a typical pattern in the other. These matters are 
the subject of ongoing study. 
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TABLE I. Results for the difference between the free en- 
ergy of hep and fee structures, as defined in Eq. (Q) with 
associated uncertainties in parentheses. Results attributed to 
Ref. Q were deduced by combining the separate results for fee 
and hep given there. PW signifies the present work. The PW 
error bounds were computed from the statistical uncertainties 
in the weights of the peaks in P(A4). 




FIG. la. Schematic representation of the close-packed 
structures. The points marked A show the positions of 
the sites in one close-packed (xy) layer; the circles show 
the boundaries of spheres occupying these sites in an ideal 
(zero-displacement) structure. The points marked B and C 
show the projections of sites in other layers (stacked along 
the z-axis) onto the xy plane; the fee and hep structures en- 
tail sequences of type ABC A . . . and ABAB . . . respectively. 
The lattice switch from fee to hep entails translations of the 
close-packed planes, as detailed in Fig. 1(b). 




FIG. lb. [Left side] The positions of the spheres in an 
arbitrary configuration of the fee structure, projected onto 
the xz plane. We show 6 layers, with 3 spheres in each; the 
sites in the top 3 layers (A,B,C) correspond to those marked 
(and underlined) A,B,C in Fig. 1(a). 

[Center] The action of the lattice switch: reading from the top, 
the first two layers (A,B) of the fee structure are invariant; 
the next two (C,A) layers are translated along the x direction 
by — 1\ and the final two (B,C) layers are translated by +t, 
where t is identified in Fig. 1(a). 

[Right side] Projections of the spheres in the resulting hep 
arrangement. Here, the displacements {«}, realizable in the 
fee structure, give two overlapping pairs of spheres (shaded) 
in the hep structure so that Ai({u}) — 2 in this case. Note 
that the picture is schematic: in particular, the density shown 
here is much lower than that chosen for the present study. 



0.03 



0.02 



0.01 



0.00 




-400 



400 



FIG. 2. The distribution of the overlap parameter M for 
a system of N=1728 spheres; the inset shows the distribu- 
tion on a logarithmic scale. The statistical uncertainties are 
smaller than the symbol size. The free-energy difference Af 
is identified from the logarithm of the ratio 1Z of the weights 
of the two peaks (Eq.(|ll|)). The smallness of Af allows both 
peaks to be displayed on one linear scale in this case. 
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